Phase-Resolved Optical Coherence Elastography: An Insight into Tissue Displacement Estimation

Robust methods to compute tissue displacements in optical coherence elastography (OCE) data are paramount, as they play a significant role in the accuracy of tissue elastic properties estimation. In this study, the accuracy of different phase estimators was evaluated on simulated OCE data, where the displacements can be accurately set, and on real data. Displacement (∆d) estimates were computed from (i) the original interferogram data (Δφori) and two phase-invariant mathematical manipulations of the interferogram: (ii) its first-order derivative (Δφd) and (iii) its integral (Δφint). We observed a dependence of the phase difference estimation accuracy on the initial depth location of the scatterer and the magnitude of the tissue displacement. However, by combining the three phase-difference estimates (Δdav), the error in phase difference estimation could be minimized. By using Δdav, the median root-mean-square error associated with displacement prediction in simulated OCE data was reduced by 85% and 70% in data with and without noise, respectively, in relation to the traditional estimate. Furthermore, a modest improvement in the minimum detectable displacement in real OCE data was also observed, particularly in data with low signal-to-noise ratios. The feasibility of using Δdav to estimate agarose phantoms’ Young’s modulus is illustrated.


Introduction
Evaluation of tissue biomechanical properties can provide important information for disease diagnosis and monitoring. Tissue stiffness has long been recognized as a biomarker of disease. Its assessment by palpation was already performed in ancient Egypt, and it is still a common tool for physical diagnosis in clinical practice [1]. Quantitative assessment of tissue biomechanics can be performed non-invasively using elastography techniques, including optical coherence elastography (OCE).
First introduced in 1998 by J. M. Schmitt [2], OCE is still an emerging technique for assessing tissue elasticity. This imaging modality combines optical coherence tomography (OCT) with a locally applied force to induce tissue displacement [3]. There are different implementations of OCE, varying the type of mechanical loading of the tissue (e.g., contact or non-contact, static, quasi-static, or dynamic, localized or global), the OCT scanning protocol, and the method used to quantify tissue deformation [3][4][5][6]. Nonetheless, they are all based on the detection of tissue deformation from consecutive OCT scans.
Speckle-tracking or phase-resolved techniques can be used to obtain tissue displacements. Speckle-tracking methods are based on intensity information and, therefore, typically require large displacements, making them unsuitable for imaging fragile tissues [7]. Phase-resolved measurements instead use phase difference information [6,7]. This provides nanoscale displacement sensitivity [3]. Combined with the microscale spatial resolution of OCT systems, OCE enables measurement scales unattainable by other elastography methods [3]. Therefore, its application in medical research has gained momentum, especially in ophthalmology [8][9][10].
The system's phase stability is a determinant for the measurement precision of phaseresolved OCE. Imperfect timing synchronization is the primary source of phase instability [11]. This problem is more relevant for swept-source systems due to the frequency jitter introduced by wavelength sweeping with mechanically moving mirrors. Several approaches have been considered and implemented to minimize the timing errors that lead to phase jumps. Commonly used solutions involve the generation of optical timing references. To correct for the nonlinearity of the swept-source frequency sweeps, Mach-Zehnder Interferometer (MZI) optical clocks are typically used. Fiber Bragg Gratings (FBG) are used as a wavelength-dependent trigger signal to synchronize the source sweep with the data digitization by the data acquisition (DAQ) board [11]. Other solutions include the implementation of a second reference arm to measure time-induced phase variations, as proposed by Vacok et al. [12]. The authors greatly improved the system's sensitivity by calibrating the phase using the recorded phase noise [12]. Recently, Li et al. described a 40-fold increase in phase stability by using a common-path OCE configuration, achieving displacement sensitivities of 0.3 nm [13]. These configurations significantly improve the phase stability of the system and, consequently, the precision of phase-resolved measurements. However, they also require changes to the instrument's setup that are not always possible or desirable.
The precision and accuracy of the method used to estimate the phase difference also play an important role in the correct assessment of tissue elasticity [6]. Several methods have been proposed to obtain the phase difference from successive interferograms. The most direct approach is obtained by subtracting the phases of consecutive interferograms acquired at the same location within a short period of time [6,7], but this is very sensitive to noise [7]. To increase the signal-to-noise ratio (SNR) and stability, a 2D autocorrelation approach has been proposed by Loupas et al. [14]. This method exploits depth-wise phase information for phase difference estimation. The size of the axial window is variable. However, to obtain an adequate SNR improvement while maintaining a high displacement axial resolution, it cannot be greater than 2% of the imaging depth [7]. Simultaneous averaging of the phase difference in the lateral and axial directions using a rectangular window has also been proposed to increase the robustness of the estimation [15]. This approach has been combined with intensity-based pixel-scale displacement tracking (speckle-tracking) to reduce additive noise and displacement-related decorrelation [15]. Phase averaging combined with a vector method has also been proposed, where incremental vectors of lateral and vertical phase variations are generated from the complex-valued signals to further reduce additive noise and eliminate the need for phase unwrapping [16,17].
The long-term goal of our research efforts is the early detection of neurodegeneration, namely the detection of Alzheimer's disease (AD) in its asymptomatic phase, by measuring changes in the mechanical properties of the retina, the visible part of the central nervous system. Significant changes in the brain's elastic properties have already been detected in AD patients [18][19][20]. Moreover, cell mechanic properties appear to be affected by the accumulation of amyloid beta that is characteristic of AD senile plaques [21]. Our goal is to detect subtle changes in retinal elasticity, highlighting the importance of accurately estimating tissue displacements. The goal is to detect subtle changes in tissue elasticity, highlighting the importance of accurately estimating tissue displacements. Reductions in Sensors 2023, 23, 3974 3 of 14 axial and temporal resolution due to the averaging of signals in depth and lateral directions must also be avoided.
In this study, we evaluated the performance of the original phase difference estimation method to better understand the influence of the magnitude of tissue displacement as well as the scatterer's depth location on the accuracy of the estimated phase difference from successive scans and to develop alternative solutions to improve the accuracy of phase difference estimation. In addition to the estimation of the phase difference from the raw interferogram data, we also evaluated the estimates from two phase-invariant mathematical manipulations of the interferogram: (i) the first-order derivative of the original interferogram and (ii) the integral of the original interferogram, and from a combination of the three approaches. The performance of each method was evaluated on numerical simulations of OCE data, with and without noise, and real data.

Theory
In phase-resolved measurements, the sample displacement ∆d is obtained by estimating the phase difference ∆ϕ [6] as: where, λ 0 is the center wavelength of the laser and n is the refractive index of the sample. As mentioned above, there are several methods to estimate ∆ϕ[z] at depth z, with different accuracies and sensitivities to noise. The A-scan intensity and phase are obtained from the measured interferogram F[k] by the inverse discrete Fourier transform (IDFT). F[k] is a discrete signal with N samples corresponding to the wavenumbers swept by the laser source from k min to k max . Hence, where j is the complex unit and the complex function f is discrete in space z. OCT A-scans are defined by the magnitude of f [z], that is, A scan [z] = | f [z]|, with z = 0, . . . , N − 1, and the phase ϕ[z] given by the angle ∠ of f [z]. Now, considering two successive interferograms, F 0 and F 1 , acquired at the same lateral location of two instants in time, t 0 and t 1 , with t 1 = t 0 + ∆t, ∆ϕ[z] can be estimated as [6,7]: where Im and Re are the imaginary and real parts of the complex function f , respectively. In this formulation, the phase difference is computed using only one evaluation of the inverse tangent function, which improves the numerical accuracy over subtracting the separately computed phases. However, it is still sensitive to noise [7]. In this work, we introduce the estimation of ∆ϕ[z] from phase-difference-invariant mathematical variations of the original interferogram: (i) the first-order derivative and (ii) the integral of the original interferogram, respectively computed using a first-order finite difference approximation of the derivative

System Setup
This study used a home-built swept-source OCE (SS-OCE) system. A schematic representation of the system is shown in Figure 1. It consists of an SS-OCT with a sweptsource laser (Axsun, Excelitas Technologies Corp., Mississauga, ON, Canada) emitting at a center wavelength of 1040 nm with a bandwidth of 110 nm and a repetition rate of 100 kHz, coupled with a piezoelectric actuator to induce tissue displacement. The output of the laser source is split by a 90 to 10 optical fiber coupler, with 90% of the light being used to generate the interferograms and the remaining 10% being directed into a FBG. The light reflected from the FBG, at a wavelength of 990 ± 1 nm, is converted to an electrical signal by a photodiode and directed to a digital delay and pulse generator (DG535, Stanford Research Systems, Sunnyvale, CA, USA) to generate trigger pulses with appropriate timing and amplitude for data acquisition. [ ], respectively, using Equation (4) after the IDFT.

System Setup
This study used a home-built swept-source OCE (SS-OCE) system. A schematic representation of the system is shown in Figure 1. It consists of an SS-OCT with a sweptsource laser (Axsun, Excelitas Technologies Corp., Mississauga, ON, Canada) emitting at a center wavelength of 1040 nm with a bandwidth of 110 nm and a repetition rate of 100 kHz, coupled with a piezoelectric actuator to induce tissue displacement. The output of the laser source is split by a 90 to 10 optical fiber coupler, with 90% of the light being used to generate the interferograms and the remaining 10% being directed into a FBG. The light reflected from the FBG, at a wavelength of 990 ± 1 nm, is converted to an electrical signal by a photodiode and directed to a digital delay and pulse generator (DG535, Stanford Research Systems, Sunnyvale, CA, USA) to generate trigger pulses with appropriate timing and amplitude for data acquisition. The light is further split 90 to 10 between the sample and reference arms to generate the interferograms. In the sample arm, the light is delivered to and collected from the sample using a 50 to 50 optical fiber coupler. An LSM03-BB objective (Thorlabs, Inc., Newton, NJ, USA), with an effective focal length of 36 mm and an entrance pupil diameter of 4 mm, was used to focus the light onto the sample. In the reference arm, the light is reflected from a fixed reference mirror. A fiber polarization controller removes polarization variations between the sample and reference signals. The resulting interference fringes, formed by coupling the sample and reference reflected light using a 50 to 50 optical fiber coupler, are detected by a balanced photodetector that subtracts the two signals to remove the common-mode noise. A DAQ board (X5-400M, Innovative Integration, Inc., Indianapolis, IN, USA) digitizes the photodetector output using the FBG-generated trigger. In addition, the inverted signal of the FBG transmission is recorded in the second analog-todigital converter (ADC) channel of the DAQ board for further post-processing jitter correction. The clock signal is provided by the internal MZI clock of the laser source.
Displacements were induced using a piezoelectric actuator (P-287, Physik Instrumente GmbH & Co. KG, Karlsruhe, Germany). The tip of the contact rod is circular, with a diameter of approximately 1.3 mm. It was positioned approximately 1.5 mm to 4.5 mm from the scan line of the OCT beam (Figure 1, inset). Mechanical vibrations were triggered by a transient pulse signal synchronized with the galvanometric mirrors and the FBGgenerated trigger to initiate data acquisition using a field-programmable gate array The light is further split 90 to 10 between the sample and reference arms to generate the interferograms. In the sample arm, the light is delivered to and collected from the sample using a 50 to 50 optical fiber coupler. An LSM03-BB objective (Thorlabs, Inc., Newton, NJ, USA), with an effective focal length of 36 mm and an entrance pupil diameter of 4 mm, was used to focus the light onto the sample. In the reference arm, the light is reflected from a fixed reference mirror. A fiber polarization controller removes polarization variations between the sample and reference signals. The resulting interference fringes, formed by coupling the sample and reference reflected light using a 50 to 50 optical fiber coupler, are detected by a balanced photodetector that subtracts the two signals to remove the common-mode noise. A DAQ board (X5-400M, Innovative Integration, Inc., Indianapolis, IN, USA) digitizes the photodetector output using the FBG-generated trigger. In addition, the inverted signal of the FBG transmission is recorded in the second analog-to-digital converter (ADC) channel of the DAQ board for further post-processing jitter correction. The clock signal is provided by the internal MZI clock of the laser source.
Displacements were induced using a piezoelectric actuator (P-287, Physik Instrumente GmbH & Co. KG, Karlsruhe, Germany). The tip of the contact rod is circular, with a diameter of approximately 1.3 mm. It was positioned approximately 1.5 mm to 4.5 mm from the scan line of the OCT beam ( Figure 1, inset). Mechanical vibrations were triggered by a transient pulse signal synchronized with the galvanometric mirrors and the FBGgenerated trigger to initiate data acquisition using a field-programmable gate array (FPGA). The width, delay, and amplitude of the transient pulse signals were controlled by an arbitrary function generator (AFG3101C, Tektronix, Beaverton, OR, USA).

Numerical Simulations
The performance of the phase difference estimators was first tested on numerical simulations of the OCE data. Based on the specifications of the swept-source laser used in our system, simulated interferograms for a single scatterer in a medium with a refractive index of 1.38 were generated using Matlab 2021a (MathWorks ® , Natick, MA, USA). Considering the spectral intensity distribution of the swept-source laser S(k) and that the refractive index n was constant with depth, the simulated interference signals I(k) were generated as [22]: where k is the wavenumber swept by the laser, ranging between k min and k max within N discrete samples, a R stands for the amplitude of the reference, a(z) is the backscatter amplitude of the object at depth z and, 2r and 2(r + z ) are the path lengths in the reference and object arms, respectively. The autocorrelation term resulting from Equation (5) was neglected. The IDFT was then used to retrieve the intensity and phase of the simulated interferogram. The scatter's depth location is precisely adjusted by changing the length of the object's arm 2(r + z ). A depth-wise shift of 90 nm in each direction was imposed to estimate the phase difference between the two interfaces. In these simulations, the interface shift was limited to ±90 nm relative to the initial location to avoid consideration of phase unwrapping.
The theoretical axial resolution of a swept-source system with the above specifications in a medium with a refractive index of 1.38 is approximately 6.7 µm. To better simulate the results obtained with this system, the initial location of the scatterer was moved 7 µm in steps of 10 nm, comprising 710 depth locations ( Figure 2). A depth range between 300 µm and 307 µm was chosen because it corresponds to the highest SNR of the system described above. For each location, the scatterer was shifted over a range of 180 nm, and the phase difference was computed. The error in phase difference estimation of each method was evaluated under ideal conditions (no noise) and after adding the system noise (real data) to the generated interferogram. The system noise was approximately 47.73 LSBrms (least significant bit root-mean-square). The original and the shifted interferograms were subjected to random noise of the same level.

Numerical Simulations
The performance of the phase difference estimators was first tested on numerical simulations of the OCE data. Based on the specifications of the swept-source laser used in our system, simulated interferograms for a single scatterer in a medium with a refractive index of 1.38 were generated using Matlab 2021a (MathWorks ® , Natick, MA, USA). Considering the spectral intensity distribution of the swept-source laser ( ) and that the refractive index was constant with depth, the simulated interference signals ( ) were generated as [22]: where is the wavenumber swept by the laser, ranging between and within discrete samples, stands for the amplitude of the reference, ( ) is the backscatter amplitude of the object at depth and, 2 and 2( + ) are the path lengths in the reference and object arms, respectively. The autocorrelation term resulting from Equation (5) was neglected. The IDFT was then used to retrieve the intensity and phase of the simulated interferogram. The scatter's depth location is precisely adjusted by changing the length of the object's arm 2( + ). A depth-wise shift of 90 nm in each direction was imposed to estimate the phase difference between the two interfaces. In these simulations, the interface shift was limited to ±90 nm relative to the initial location to avoid consideration of phase unwrapping.
The theoretical axial resolution of a swept-source system with the above specifications in a medium with a refractive index of 1.38 is approximately 6.7 μm. To better simulate the results obtained with this system, the initial location of the scatterer was moved 7 μm in steps of 10 nm, comprising 710 depth locations ( Figure 2). A depth range between 300 μm and 307 μm was chosen because it corresponds to the highest SNR of the system described above. For each location, the scatterer was shifted over a range of 180 nm, and the phase difference was computed. The error in phase difference estimation of each method was evaluated under ideal conditions (no noise) and after adding the system noise (real data) to the generated interferogram. The system noise was approximately 47.73 LSBrms (least significant bit root-mean-square). The original and the shifted interferograms were subjected to random noise of the same level.

Phase Difference Measurements
All methods were implemented using Matlab 2021a (MathWorks ® , Natick, MA, USA). The workflow is shown in Figure 3. The original phase difference estimate (Δ ) was obtained after the IDFT of each interferogram to obtain the respective A-scan magnitude and phase. The phase difference between consecutive A-scans was then calculated using Equation (4) ( Figure 3A). For the proposed phase difference estimators, the first-

Phase Difference Measurements
All methods were implemented using Matlab 2021a (MathWorks ® , Natick, MA, USA). The workflow is shown in Figure 3. The original phase difference estimate (∆ϕ ori ) was obtained after the IDFT of each interferogram to obtain the respective A-scan magnitude and phase. The phase difference between consecutive A-scans was then calculated using Equation (4) ( Figure 3A). For the proposed phase difference estimators, the first-order derivative (F d [k]) and integral (F int [k]) of the original interferogram were computed prior to the IDFT. An approximation of F d [k] was obtained using the difference between the discrete channels ( Figure 3B). For an approximation of F int [k], the cumulative integral using the trapezoidal method was used ( Figure 3C). Following the IDFT, ∆ϕ d and ∆ϕ int were estimated from Equation (4). Phase unwrapping was applied to the estimated phase differences ∆ϕ ori , ∆ϕ d , and ∆ϕ int , affecting the displacement estimation error of all methods equally. Equation (1) was used to calculate the displacements ∆d ori , ∆d d , and ∆d int , from the estimated phase differences ∆ϕ ori , ∆ϕ d , and ∆ϕ int , respectively. between the discrete channels ( Figure 3B). For an approximation of [ ], the cumulative integral using the trapezoidal method was used ( Figure 3C). Following the IDFT, Δ and Δ were estimated from Equation (4). Phase unwrapping was applied to the estimated phase differences Δ , Δ , and Δ , affecting the displacement estimation error of all methods equally. Equation (1) was used to calculate the displacements Δ , Δ , and Δ , from the estimated phase differences Δ , Δ , and Δ , respectively.

Static Conditions
The OCT interference signal of a highly reflective sample (gold-coated mirror) was used to determine the temporal stability of the system and the performance of each estimator under static (real) conditions, i.e., without inducing any displacement. The influence of the imaging depth, and consequently the SNR, was evaluated by changing the relative depth location of the sample. Data were recorded at optical path differences (OPDs) between 0 µm and 3000 µm in steps of 500 µm, with OPD = 0 µm being the depth location with the highest measured SNR (71.19 dB) and OPD = 3000 μm the lowest (40.16 dB).
The minimum detectable displacements were calculated to determine the resolution of each estimator. Fundamental limitations of the minimum detectable displacement in phase-sensitive measurements are phase stability and the precision of the Δ computation. The former is limited by the system's temporal stability and is inherently related to the instrument's setup. The measured temporal stability of the system was 396.9 ± 46.7 [23]. No changes were observed in the system's temporal stability with a decrease in the system SNR.
The displacement resolution is linearly dependent on the standard deviation of the measured phase difference ( ), as shown in Equation (1) [6]. Thus, the was calculated as a metric of the minimum detectable displacement. The lower the error associated with the displacement estimation, the higher the accuracy. The root-mean-square error (RMSE) was computed to provide additional information about the accuracy of the displacement estimation.

Dynamic Conditions
Data were collected using the Motion-Brightness (M-B) scanning protocol [3]. Axial scans were repeated 768 times at a fixed spatial location to track surface motion at high  (4):

Static Conditions
The OCT interference signal of a highly reflective sample (gold-coated mirror) was used to determine the temporal stability of the system and the performance of each estimator under static (real) conditions, i.e., without inducing any displacement. The influence of the imaging depth, and consequently the SNR, was evaluated by changing the relative depth location of the sample. Data were recorded at optical path differences (OPDs) between 0 µm and 3000 µm in steps of 500 µm, with OPD = 0 µm being the depth location with the highest measured SNR (71.19 dB) and OPD = 3000 µm the lowest (40.16 dB).
The minimum detectable displacements were calculated to determine the resolution of each estimator. Fundamental limitations of the minimum detectable displacement in phase-sensitive measurements are phase stability and the precision of the ∆ϕ computation. The former is limited by the system's temporal stability and is inherently related to the instrument's setup. The measured temporal stability of the system was 396.9 ± 46.7 ps [23]. No changes were observed in the system's temporal stability with a decrease in the system SNR.
The displacement resolution is linearly dependent on the standard deviation of the measured phase difference (σ ∆ϕ ), as shown in Equation (1) [6]. Thus, the σ ∆ϕ was calculated as a metric of the minimum detectable displacement. The lower the error associated with the displacement estimation, the higher the accuracy. The root-mean-square error (RMSE) was computed to provide additional information about the accuracy of the displacement estimation.

Dynamic Conditions
Data were collected using the Motion-Brightness (M-B) scanning protocol [3]. Axial scans were repeated 768 times at a fixed spatial location to track surface motion at high frame rates, corresponding to a scan time of 7.68 ms per location. B-scans were generated by temporally aligning the different spatial locations. Imaging was performed at 256 lateral locations. Data were acquired over an approximately 3 mm line in steps of 11.72 µm (Figure 1, inset).
Transient pulses were used to induce motion in homogeneous agarose phantoms. The surface motion was generated once at each spatial location, with a delay of 100 µs (corresponding to 10 interferograms), using transient pulses with widths between 100 µs and 600 µs, in steps of 100 µs, with a fixed amplitude of 50 mV. A total of three measurements were made for each pulse width. The average laser power at the sample arm was approximately 1.06 mW.
The estimator with the highest accuracy in both simulated and static data was then used to obtain the phase difference at the phantom boundary. The surface displacement was calculated using Equation (1) to derive the elastic modulus (Young's modulus) of the phantom, E = 3ρC 2 s , where ρ is the material density and C s is the shear wave velocity. The refractive index of agarose phantoms was calculated from their concentration c as n p = 0.0014c + 1.333 [24].
The propagation velocity of the generated surface waves (Rayleigh waves) C R was obtained from the propagation distance of a transient pulse ∆D during the time ∆t as C R = ∆D/∆t. The C s was then obtained using the correlation C s = 1.05C R [7]. Average cross-correlations between tissue displacement curves over time (n = 10) in 0.5 mm increments were used to obtain ∆D.

Phantom Preparation
Homogeneous agarose phantoms, with tissue-mimicking properties, were prepared by gently mixing agar with distilled water at a high temperature (85 • C) in concentrations of 10 g/L, 15 g/L, and 25 g/L while stirring. Glass microspheres, with diameters of 152 ± 32 µm, were added to the mixture at a concentration of 10 g/L to increase optical scattering. The mixture was then placed in containers for molding and cured for approximately 24 h. The resulting phantoms were cylindrical, with diameters and heights of approximately 5.2 cm and 3.5 cm, respectively.

Statistical Analysis
Statistical differences between the squared errors of the ∆ϕ estimators were assessed using the nonparametric Wilcoxon test for paired samples after the Kolmogorov-Smirnov test was used to reject the normality hypothesis. Differences were considered significant at a significance level of 0.05.

Simulation Results
A simulation of the interferogram expected from our OCE system was generated to better calculate the error associated with phase difference estimation using different methods. The initial location of the scatterer varied over 7 µm (710 depth locations).
Each method's performance depended on the scatterer's initial depth location. The error, measured as the difference between simulated and measured displacements associated with ∆d ori , ∆d d , and ∆d int , was lower for a total of 267, 204, and 239 depth locations, respectively, out of the 710 depth locations tested. This shows that, depending on the initial depth location of the scatterer, ∆d ori will be a better approach 37.61% of the time, while ∆d d will be a better approach 28.78% of the time, and ∆d int will be a better approach the remaining 33.66% of the time. Figure 4 shows examples of the error associated with displacement estimation at three different depths, where either ∆d ori (Figure 4A), ∆d d ( Figure 4B), or ∆d int ( Figure 4C) estimators yielded values closer to the real ones (lower error).  As expected, when zero displacements (static) were simulated, the error associated with the displacement calculation was zero for all methods. Increasing the magnitude of the simulated displacement in any direction decreased the accuracy of all methods. Interestingly, although the error associated with displacement computation increased with distance for all methods, they showed different trends. As shown in Figure 4, the estimation errors for Δ and Δ and Δ had opposite and complementary trends. As the magnitude of the displacement induced increased, Δ tended to underestimate the displacement, while Δ and Δ tended to overestimate it. In view of this, we propose a novel displacement estimator Δ computed from Δ which combines all the phase difference estimation approaches as Δ = (Δ + Δ + Δ ) 3 ⁄ . The error associated with the estimation of Δ is also shown in Figure 4. As observed, Δ provides a more accurate prediction of the displacement than either approach alone.
The RMSE was used to evaluate the quality of each estimator. The kernel density distributions of the RMSE of the displacement computation for all depth locations, obtained with Δ and Δ in simulations without and with noise, are shown in Figure  5. The proposed approach performs better, with a pronounced decrease in the median error and interquartile range (IQR). Detailed values are given in Table 1.  As expected, when zero displacements (static) were simulated, the error associated with the displacement calculation was zero for all methods. Increasing the magnitude of the simulated displacement in any direction decreased the accuracy of all methods. Interestingly, although the error associated with displacement computation increased with distance for all methods, they showed different trends. As shown in Figure 4, the estimation errors for ∆d ori and ∆d d and ∆d int had opposite and complementary trends. As the magnitude of the displacement induced increased, ∆d ori tended to underestimate the displacement, while ∆d d and ∆d int tended to overestimate it. In view of this, we propose a novel displacement estimator ∆d av computed from ∆ϕ av which combines all the phase difference estimation approaches as ∆ϕ av = (∆ϕ ori + ∆ϕ d + ∆ϕ int )/3. The error associated with the estimation of ∆d av is also shown in Figure 4. As observed, ∆d av provides a more accurate prediction of the displacement than either approach alone.
The RMSE was used to evaluate the quality of each estimator. The kernel density distributions of the RMSE of the displacement computation for all depth locations, obtained with ∆d ori and ∆d av in simulations without and with noise, are shown in Figure 5. The proposed approach performs better, with a pronounced decrease in the median error and interquartile range (IQR). Detailed values are given in Table 1. As expected, when zero displacements (static) were simulated, the error associated with the displacement calculation was zero for all methods. Increasing the magnitude of the simulated displacement in any direction decreased the accuracy of all methods. Interestingly, although the error associated with displacement computation increased with distance for all methods, they showed different trends. As shown in Figure 4, the estimation errors for Δ and Δ and Δ had opposite and complementary trends. As the magnitude of the displacement induced increased, Δ tended to underestimate the displacement, while Δ and Δ tended to overestimate it. In view of this, we propose a novel displacement estimator Δ computed from Δ which combines all the phase difference estimation approaches as Δ = (Δ + Δ + Δ ) 3 ⁄ . The error associated with the estimation of Δ is also shown in Figure 4. As observed, Δ provides a more accurate prediction of the displacement than either approach alone.
The RMSE was used to evaluate the quality of each estimator. The kernel density distributions of the RMSE of the displacement computation for all depth locations, obtained with Δ and Δ in simulations without and with noise, are shown in Figure  5. The proposed approach performs better, with a pronounced decrease in the median error and interquartile range (IQR). Detailed values are given in Table 1.   As expected, a decrease in performance was observed for both methods when noise was added to the generated interferogram. Nevertheless, the proposed estimator (∆d av ) performs significantly better than the original one (∆d ori ).

Minimum Detectable Displacement
Both ∆d ori and ∆d av were used to estimate the displacement from the same elastography data in static conditions (real data) at different OPDs and their corresponding minimum detectable displacements. The kernel density distributions of the RMSE associated with the displacement estimation at 256 different lateral locations using ∆d ori and ∆d av are shown in Figure 6A. Figure 6B shows the normalized RMSE of all locations for all OPDs. As demonstrated by the numerical simulation, the error associated with displacement estimation was equivalent for ∆d ori and ∆d av for small displacements, and both converged to zero when no displacement was induced ( Figure 4). Thus, for static measurements, equivalent error distributions were expected for both methods ( Figure 6A). Nevertheless, a significant decrease in the overall RMSE was obtained for OPDs larger than 1000 µm (lower SNR). performs significantly better than the original one (Δ ). Table 1. Detailed values of the displacement root-mean-square error (RMSE) distribution in simulated interferograms without and with noise using the original (Δ ) and the proposed estimator (Δ ). Data with noise are presented as the mean ± standard deviation. * Statistically significant at p < 0.001.

Minimum Detectable Displacement
Both Δ and Δ were used to estimate the displacement from the same elastography data in static conditions (real data) at different OPDs and their corresponding minimum detectable displacements. The kernel density distributions of the RMSE associated with the displacement estimation at 256 different lateral locations using Δ and Δ are shown in Figure 6A. Figure 6B shows the normalized RMSE of all locations for all OPDs. As demonstrated by the numerical simulation, the error associated with displacement estimation was equivalent for Δ and Δ for small displacements, and both converged to zero when no displacement was induced ( Figure 4). Thus, for static measurements, equivalent error distributions were expected for both methods ( Figure  6A). Nevertheless, a significant decrease in the overall RMSE was obtained for OPDs larger than 1000 μm (lower SNR).  Table 2). For OPDs equal to and greater than 1000 μm, an improvement in the accuracy of the displacement estimation was Figure 6. Kernel density estimates of the root-mean-square error (RMSE) displacement distribution for 256 depth locations (A) and mean normalized RMSE (B) of a static gold mirror at multiple optical path differences (OPD) using the original (∆d ori ) and the proposed estimator (∆d av ). The median (dashed) and the first and third quartiles (dotted) are shown. * Statistically significant at p < 0.001. The minimum detected displacement at different OPDs was obtained from σ ∆ϕ ori and σ ∆ϕ av at 256 different locations using Equation (1) ( Table 2). For OPDs equal to and greater than 1000 µm, an improvement in the accuracy of the displacement estimation was achieved using ∆d av estimation. Improvements between 0.01% and 0.37% were observed depending on the imaging depth. No significant changes were observed for optimal SNR (OPD = 0 µm).

Phantom Displacement
The proposed displacement estimator (∆d av ) was used to infer the elastic properties of agarose phantoms with different elastic properties. Figure 7 shows the spatiotemporal displacement map for a phantom with an agarose concentration of 10 g/L after a transient pulse of 500 µs ( Figure 7A), displacement curves at three lateral locations at progressively higher distances from the piezoelectric actuator (0, 1.5, and 3 mm) ( Figure 7B), and the structural B-scan of the phantom's interface superimposed with ∆d av at three times (3, 4, and 5 ms).

Phantom Displacement
The proposed displacement estimator (Δ ) was used to infer the elastic properties of agarose phantoms with different elastic properties. Figure 7 shows the spatiotemporal displacement map for a phantom with an agarose concentration of 10 g/L after a transient pulse of 500 μs ( Figure 7A), displacement curves at three lateral locations at progressively higher distances from the piezoelectric actuator (0, 1.5, and 3 mm) ( Figure 7B), and the structural B-scan of the phantom's interface superimposed with Δ at three times (3, 4, and 5 ms).  The propagation of the Rayleigh wave velocity was calculated using the displacement estimation based on the ratio between the propagation distance and the time to reach the maximum displacement. No discernible differences were observed when using transient pulses of different widths. Therefore, all 18 measurements were used to obtain the shear wave propagation velocity of each phantom. As expected, a linear increase in shear wave velocity was observed with an increase in the phantom's agarose concentration ( Figure 8A), which corresponded to an increase in the phantom's Young's modulus ( Figure 8B).
Sensors 2023, 23, 3974 11 of The propagation of the Rayleigh wave velocity was calculated using the displaceme estimation based on the ratio between the propagation distance and the time to reach t maximum displacement. No discernible differences were observed when using transie pulses of different widths. Therefore, all 18 measurements were used to obtain the she wave propagation velocity of each phantom. As expected, a linear increase in shear wa velocity was observed with an increase in the phantom's agarose concentration (Figu 8A), which corresponded to an increase in the phantom's Young's modulus (Figure 8B

Discussion
Quantitative evaluation of tissue elastic properties by phase-resolved OCE depen on the measurement noise, the precision of the phase difference estimation and hence t tissue displacement, and the method used to derive the tissue's Young's modulus. Me urement noise is intrinsically related to the phase stability of the system. This can be i proved by modifying the instrument setup, typically by introducing timing references. this study, an MZI clock and FBG-triggering were used to improve timing synchroni tion and achieve temporal stabilities consistent with those reported in the literatu [11,13,25]. Various approaches have also been proposed to reduce noise sensitivity a improve the accuracy of phase difference estimation. However, these rely heavily depth (often combined with lateral) averaging of the interferogram data, which redu their resolution.
In this study, we characterize the performance of phase difference estimation a function of the magnitude of the tissue displacement and the depth location of the sc terer. The accuracy of the original phase difference estimation method (Δ ) was eva ated, as well as that of phase-difference estimators derived from: (i) the first-order der ative (Δ ), (ii) the integral (Δ ) of the interferogram data, and (iii) the combination three phase-difference estimators (Δ ). These estimators were tested on simulated O data in which the displacements could be precisely adjusted and were known a priori. W observed that the initial depth location of the scatterer affected the precision of the Δ Δ , and Δ estimators. Each approach was the most accurate in 1/3 of the initial dep locations tested. Since the initial location of the scatterer is not known from real measu ments, we predicted that the best displacement estimation would be achieved using combination of all approaches.
We demonstrated that Δ was the most accurate displacement estimator in t simulated data. For small displacements, both methods achieved similar accuracy. Ho ever, as the simulated displacement increased, the error associated with the estimati was lower for Δ than for Δ . This was possible due to the complementary and o posite trends of the Δ estimation and those from phase-invariant mathematical va ations of the original interferogram (Δ and Δ ). While the latter overestimated t

Discussion
Quantitative evaluation of tissue elastic properties by phase-resolved OCE depends on the measurement noise, the precision of the phase difference estimation and hence the tissue displacement, and the method used to derive the tissue's Young's modulus. Measurement noise is intrinsically related to the phase stability of the system. This can be improved by modifying the instrument setup, typically by introducing timing references. In this study, an MZI clock and FBG-triggering were used to improve timing synchronization and achieve temporal stabilities consistent with those reported in the literature [11,13,25]. Various approaches have also been proposed to reduce noise sensitivity and improve the accuracy of phase difference estimation. However, these rely heavily on depth (often combined with lateral) averaging of the interferogram data, which reduces their resolution.
In this study, we characterize the performance of phase difference estimation as a function of the magnitude of the tissue displacement and the depth location of the scatterer. The accuracy of the original phase difference estimation method (∆ϕ ori ) was evaluated, as well as that of phase-difference estimators derived from: (i) the first-order derivative (∆ϕ d ), (ii) the integral (∆ϕ int ) of the interferogram data, and (iii) the combination of three phase-difference estimators (∆d av ). These estimators were tested on simulated OCE data in which the displacements could be precisely adjusted and were known a priori. We observed that the initial depth location of the scatterer affected the precision of the ∆d ori , ∆d d , and ∆d int estimators. Each approach was the most accurate in 1/3 of the initial depth locations tested. Since the initial location of the scatterer is not known from real measurements, we predicted that the best displacement estimation would be achieved using a combination of all approaches.
We demonstrated that ∆d av was the most accurate displacement estimator in the simulated data. For small displacements, both methods achieved similar accuracy. However, as the simulated displacement increased, the error associated with the estimation was lower for ∆d av than for ∆d ori . This was possible due to the complementary and opposite trends of the ∆d ori estimation and those from phase-invariant mathematical variations of the original interferogram (∆d d and ∆d int ). While the latter overestimated the simulated displacement with increasing displacement, the former underestimated it. Overall, a reduction in the median RMSE associated with displacement prediction of approximately 70% and 85% was observed for simulated data without and with noise, respectively. Despite the relatively modest gain, our goal is to use OCE to detect changes in the mechanical properties of the retina that may be associated with Alzheimer's disease and other neurodegenerative disorders before structural changes can be detected. Regardless of its magnitude, any correction will bring us closer to our goal.
No significant changes were expected for static OCE data, as both methods converge to zero in the absence of displacement. However, ∆d av resulted in an improvement in the displacement estimation depending on the initial relative depth location of the sample and, consequently, the SNR. The ∆d av resulted in a lower RMSE than ∆d ori for imaging depths at OPDs greater than 1000 µm, where the SNR is lower. Comparatively lower minimum detectable displacements could be measured in signals with low SNR. Thus, we have shown that the proposed method can significantly improve the resolution of displacement estimation from phase data. It is noteworthy that the typical correlation between the measured phase difference and the OCT SNR (σ ∆ϕ ≈ 1/ √ SNR) [6] was not found. However, since we aim to compare the accuracy of both methods, the absolute values of the minimum detectable displacements do not have a direct impact.
The feasibility of the proposed method to estimate tissue displacement under dynamic conditions was also evaluated. Based on ∆d av , the shear wave velocity was estimated for phantoms with different elastic properties. As expected, an increase in agarose concentration resulted in an increase in the shear wave velocity and Young's modulus. Although the benefits of increased displacement measurement accuracy in wave-based OCE may not be readily apparent, they can have a significant impact on elasticity estimation when considering numerical methods. An example is the application of Newton's method to a wave solution represented by the method of fundamental solutions [26]. We have recently shown that, when using Newton's-type iterations to recover the elastic properties of the medium, 1% noise in the data can lead to a two-to three-orders-of-magnitude higher error in the recovered elastic coefficients [26].
We have illustrated that we can develop novel solutions to improve the accuracy of phase estimation by considering external factors such as the magnitude of the tissue displacement and the depth location of the scatterer. We have also illustrated that the proposed estimator, obtained by combining three estimators for the phase difference, although modest in gain, improves the accuracy of displacement estimation without the loss of temporal or axial resolution. The main drawback of the proposed method is the increase in computational time. However, it is important to mention that the improvement comes at no cost since no changes to the optical setup or the instrument are required. The proposed method is also easy to implement, independent of the system setup, and can even be applied to previously recorded data.  Data Availability Statement: The code and raw data supporting the conclusions of this manuscript will be made available by the authors upon formal and reasonable request.